Spatial Analysis

Load data

Maps data

We use three different shapefiles for the continental U.S. land mass, the State waters of maine, new hampshire, massachusets, connecticut, rhode island, new york, new jersey, delaware, maryland, virginia, north carolina and the North East EEZ.

1.- Land shapefile; covers the US land territory for visualization. Data provided from the map package.

2.- State waters; covers the state waters of the NE US states. Data from data.gov.

License: No license information was provided. If this work was prepared by an officer or employee of the United States government as part of that person’s official duties it is considered a U.S. Government Work.

3.- EEZ shapefile; Used the Sea Around Us shapefle updated to June 2016.

unique(land_sf$ID)

Spatial component

Create base shapefile for computation

As a first step we need to divide the NE US EEZ among the different states. For that we expanded state waters up to the 200 nautical miles to then estimate the percentage that each expanded-state-waters occupied. Note that in all cases these areas overlapped and percentages accounted for it. We did this by following these steps:

• 1. Expand state waters using a buffer
• 2. Grid that buffer on a 0.3 resolution
• 3. Crop the buffer to the EEZ

1. Expand Spatial regions (a.k.a buffers)

We set a buffer of 410000 m (410 km, ~ 221 nm) that overshoots the EEZ a bit, but is eventually cropped

# Buffer state waters
state_bf = st_buffer(state_sf, 410000) %>% 
  st_transform(4326) %>% # to match shape
  st_set_crs(4326)

# ------------------------------ #
# Testing and visualizing the buffer 
# ------------------------------ #

state_bf_plot <- ggplot() +
  geom_sf(data = eez_sf, aes(), fill = NA) +
  geom_sf(data = state_bf, aes(color = state), fill = NA)
#   

# ggsave("../Results/Partial/state_waters_buffer.png",state_bf_plot)

# ------------------------------ #

state_bf_plot

2. Expand grid within buffer

Here we expand a grid within the buffer so we can estimate the proportion of each state

2.1 Merge grid and buffers

Once we have a gridded area, we converted the grid to a sf so we can merge it with the buffered states and finally filter out everything outside the states polygon

# ---------------- #
# Convert grid to sf
# ---------------- #
grid_sf <- st_as_sf(ne_grid,
             coords = c("lon", "lat"),
             crs = 4326) %>% 
  st_join(state_bf) %>% 
  filter(!is.na(state))

# Create data frame for future computations
# Note, will be used in next chunk
grid_bf_dt <- as.data.frame(grid_sf) %>%
    select(index,state)
    # group_by(state) %>% 
    # summarise(length(index))


# ---------------- #
# [TEST] 
# Visualization of grid
# ---------------- #

gridExtra::grid.arrange(
  # Overall (overlapping) position
  ggplot(grid_sf) +
    geom_sf(aes(color = state), alpha = 0.3) +
    geom_sf(data = eez_sf, aes(),fill =NA) + 
    theme(legend.position = ""),
  # Showing each state separatley
  ggplot() +
    geom_sf(data = grid_sf, aes(color = state),size = 0.1, alpha = 0.3) +
    facet_wrap(~state) +
    theme(legend.position = "top"),
  nrow = 1)

3. Crop buffers to EEZ

Finally, we crop the grided buffers to within the EEZ to capture the actual water space.

Note: This step takes quite a while because of the size of the EEZ shapefile. No, you cannot use st_simplify() here

Interpolation rutine

Here we interpolate the survey data within the previously created grid.

• We removed cases where wtcpue < 0

Function for interpolation

This is the main function used to interpolate the data per year. It follows a Triangular Irregular Surface method using the interp::interp() function.

tis <- function(input_data, grid, yr, taxa, reg){
  
  # --------------- #
  # Testing
  # print(paste(yr))
  # yr = 1976
  # --------------- #
  
  # Filter data
  data <- input_data %>% 
    filter(year == yr,
           spp == taxa,
           region == reg
    ) %>% 
    # Average duplicated hauls in the same spot
    group_by(region,year,spp,lat,lon) %>%
    summarise_at(vars(wtcpue),
                 mean,na.rm = T) %>%
    filter(wtcpue > 0)
    
  # Only interpolate cases where there is more than 3 rows
  # Marked by the function 
    if(nrow(data) <= 3){
      fit_tin <- tibble()
    }else{
      
      # Triangular Irregular Surface
      fit_tin <- interp::interp( # using {interp}
        x = data$lon,           # the function actually accepts coordinate vectors
        y = data$lat,
        z = data$wtcpue,
        xo = grid$lon,     # here we already define the target grid
        yo = grid$lat,
        output = "points"
      ) %>% 
        bind_cols() %>% 
        bind_cols(grid) %>%
        mutate(year = yr,
               region = reg,
               spp = taxa) %>% 
        select(index, state, year, region, spp, lon=x, lat=y, value = z)
      
    }
   
    return(fit_tin)
}

# Test me
# Needs variables in Control panel
# Test no data: "Alosa aestivalis", reg = "Northeast US Fall", yr = 1974 
# tis(input_data = ocean_data, grid = grid_eez_df, taxa = "Illex illecebrosus", reg = "Northeast US Fall", yr = 1973)



# lapply(years,tis,input_data = ocean_data, grid = grid_eez_df, taxa = "Illex illecebrosus", reg = regions[2])

Run function

Ocean adapt data

• Using only the Northeast US Fall and Spring bottom trawl survey data for now

Control Pannel

Load required data. Note that some of it has been previously created

head(spp)

[1] "Alosa aestivalis"     "Alosa pseudoharengus"
[3] "Alosa sapidissima"    "Amblyraja radiata"   
[5] "Ammodytes dubius"     "Anarhichas lupus"    

Run routine

# single species run
# run_tis(input_data = ocean_data,
#         grid = grid_eez_df,
#         years = years,
#         reg = regions,
#         taxa = "Illex illecebrosus"
# )


# Run them all in parallel
parallel::mclapply(spp,
                   run_tis, 
                   input_data = ocean_data,
                   grid = grid_eez_df,
                   years = years,
                   reg = regions
)

Results

Points

Map

Area plot

This graph shows the proportion of each State over time

Area plot (running)

This graph shows the proportion of each State smoothed over a 10 years running mean. It helps seeing the trends better (I Think…)